Influence of Island Diffusion on Submonolayer Epitaxial Growth 
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We investigate the kinetics of submonolayer epitaxial growth which is driven by a fixed flux 
of monomers onto a substrate. Adatoms diffuse on the surface, leading to irreversible aggregation 
of islands. We also account for the effective diffusion of islands, which originates from hopping 
processes of their constituent adatoms, on the kinetics. When the diffusivity of an island of mass k 
scales as the (mean-field) Smoluchowski rate equations predicts steady behavior for < fi < 1, 
with the concentration Ck of islands of mass k varying as fc -1 - 3- ^/ 2 . For fj, > 1, a quasi-static 
approximation to the rate equations predicts slow continuous evolution in which the island density 
increases as (lnt) M//2 . A more refined matched asymptotic expansion reveals unusual multiple-scale 
mass dependence for the island size distribution. Our theory also describes basic features of epitaxial 
growth in a more faithful model of growing circular islands. For epitaxial growth in an initial 
population of monomers and no external flux, a scaling approach predicts power-law island growth 
and a mass distribution with a behavior distinct from that of the non-zero flux system. Finally, 
we extend our results to one- and two-dimensional substrates. The physically-relevant latter case 
exhibits only logarithmic corrections compared to the mean-field predictions. 

PACS numbers: 68.35.Fx, 36.40.Sx, 66.30.Fq, 82.20.Wt, 05.40.+j 



I. INTRODUCTION 

Submonolayer epitaxial thin film growth involves de- 
position of atoms onto a substrate and diffusion of these 
adatoms (monomers) leading to their aggregation into 
islands of ever-increasing size (!],§]. The resulting island 
morphology and mass distribution ultimately depend on 
these diffusion processes on the substrate. While it has 
long been recognized that these mass transport details 
are crucial to epitaxial growth Q, its ramifications are 
still incompletely understood. Part of the reason for slow 
progress is that a variety of microscopic details can and 
do influence the rate at which adsorbate islands move 

One direction of previous investigation was based on a 
picture that only monomers can diffuse while larger is- 
lands are immobile |^ In this case, the island density 
N scales with time as T 1 / 3 before it reaches a maximum 
density which scales with flux F as F 1 / 3 . Subsequent 
work revealed, however, that basic results are sensitive 
to minor alterations to the mass transport mechanism. 
For example, if both monomers and dimers diffuse while 
larger islands remain immobile, the exponents of the time 
and flux dependences change M. More generally, these 
exponents depend on the threshold size between mobile 
and immobile islands. 

On the other hand, there has been increasing apprecia- 
tion that adatom hopping continues to occur even 
when they are incorporated into islands of arbitrary size. 
This leads to a non-zero diffusivity of such islands. For 
appropriate systems and experimental conditions there 
is ample evidence that the effective island diffusivity 
has a power-law dependence on mass k 10 - 13 1 . This 
adatom hopping also typically leads to islands maintain- 
ing compact shapes, an important simplifying feature for 



theoretical modeling. 

The goal of this paper is to determine the aggregation 
kinetics of compact islands in the submonolayer regime 
with power law and more general mass-dependent island 
diffusivities [p"5| . An essential, and at first sight surpris- 
ing feature about such compact islands is that their re- 
activity depends only logarithmically on island size in 
two dimensions. Accordingly, our theoretical treatment 
is based on a model of point-like islands. Such a descrip- 
tion should accurately describe island growth, up to these 
logarithmic corrections in the low coverage limit. 

Within the mean-field Smoluchowski rate equations, 
we will consider a system with: (i) fixed monomer flux, 
(ii) point-like islands with diffusivity which decays faster 
than inversely with island mass, and (iii) irreversible 
mass-conserving coalescence of islands. We will show 
that the number of islands grows logarithmically in time 
and that the island mass distribution exhibits a multiple- 
scale mass dependence. Our predictions are also found 
to apply to a more faithful model of epitaxial growth in 
which islands are growing circles, with the radius of a 
mass-fc island proportional to y/k. Thus the point-island 
model provides a useful framework to describe submono- 
layer epitaxial growth and gives quantitative predictions 
which are robust to variations in model parameters. 

In Sec. II, we introduce our model and discuss the 
applicability of the Smoluchowski approach to epitaxial 
growth. In Sec. Ill, we present our main results about the 
asymptotic growth of islands based on the rate equations. 
Different behaviors arise depending on whether the mo- 
bility exponent fi is smaller, larger, or equal to unity. We 
also discuss the flux dependence of the maximum island 
density and the time range where our theory should ap- 
ply. In Sec. IV, we briefly discuss island growth kinetics 
in the absence of external flux. In Sec. V, we present nu- 
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merical simulations and compare results for the epitaxial 
growth of point islands and growing circular islands with 
radius proportional to \fk. Sec. VI contains a summary 
and discussion. In Appendix A, we generalize our the- 
ory to one-dimensional substrates and in Appendix B we 
outline a more accurate treatment for two-dimensional 
substrates which accounts for the logarithmic corrections 
in the reaction rate. 



II. MICROSCOPIC MODEL 

We consider submonolayer epitaxial growth which re- 
sults from the irreversible deposition of atoms from a gas 
onto a substrate and subsequent irreversible aggregation. 
Adatoms (or monomers) diffuse on the substrate and 
aggregate upon colliding, creating dimers and larger is- 
lands. We consider the class of models for which adatoms 
which are already incorporated into islands can continue 
to diffuse - this often occurs along the periphery of is- 
lands (see e.g., Ref. fHsfl). Such adatom hopping induces 
an effective diffusion in which islands of mass k hop on 
the substrate with diffusion coefficient D k . This micro- 
scopic mechanism also typically leads to compact island 
shapes. Whenever two islands meet and aggregate, we 
assume that adatom hopping quickly causes the result- 
ing aggregate to become compact. We therefore treat 
islands as circular throughout our modeling. 

Within this picture for islands, a classic way to calcu- 
late the island mass distribution is based on the Smolu- 
chowski rate equations [|l6Ul7| . This approach requires 
knowledge of the rate at which an island of mass i and 
an island of mass j aggregate to form an island of mass 
i+j. If the aggregation process is diffusion-controlled and 
islands are spherical, this aggregation rate is given by the 
Smoluchowski formula ~ (-D, + Dj)(Ri + Rj) d ~ 2 in 
d dimensions . Here R j is the radius of an island of 
mass i. This formula applies for d > 2, while in the physi- 
cally relevant case of two-dimensional substrates, there is 
only a slow logarithmic dependence on island radius [ fL7| . 
Thus a reasonable starting point for theoretical investi- 
gation is to ignore the logarithmic term; this significantly 
simplifies the resulting analysis. In Appendix B we will 
return to the two-dimensional case and show that these 
logarithmic corrections do not alter our main findings but 
rather give rise to a logarithmic rcnormalization of the 
monomer flux. 

As discussed in the Introduction, we consider the is- 
land diffusion coefficient to be a homogeneous function 
of island mass, Dk ~ kr^ with fj, non-negative on basic 
physical grounds. In appropriate time units the reaction 
rate is 

^=r"+r". (i) 

However, our approach can equally-well be applied to any 
functional dependence of Dk on k which decays faster 
than k~ x as k — > oo. 



For a self-contained discussion, we give a qualitative 
argument for the exponent value fj, = 3/2 for the "pe- 
riphery" adatom hopping mechanism; this has been ob- 
tained previously by a Langevin approach |]l3| , |l8| . In pe- 
riphery diffusion, adatoms on the edge of an island can 
hop freely to neighboring sites on the periphery. Con- 
sider an island of linear size R 3> 1. In a time interval 
At ~ R 2 an adatom on the edge will typically explore the 
entire island boundary, so this adatom will typically be 
a distance R from its initial location. Hence the effective 
island center-of-mass displacement is 6 x ~ R/R 2 = R^ 1 . 
If each edge adatom diffuses independently, the total 
center-of-mass displacement Ax is the sum of R indepen- 
dent identically distributed random variables with vari- 
ance Sx ~ R~ l . This implies Ax ~ */ R(8x) 2 ~ R~ x l 2 . 
Thus the effective center-of-mass diffusion coefficient is 
D R ~ (Ax) 2 /At ~ R- 3 , or D k ~ fc~ 3/2 . For a d- 
dimensional substrate, a straightforward generalization 
of this argument gives fj, = 1 + l/d. Similarly, in the 
case of so-called "terrace" diffusion [ fL3"| , it is found that 
Dk ~ k^ 1 independent of d. 

Before presenting detailed results, we discuss some lim- 
itations of our rate equation approach. As stated previ- 
ously, we ignore the effect of a finite island radius on the 
form of the reaction rates, as this dependence is only loga- 
rithmic in two dimensions. Moreover, as the coverage in- 
creases a rate equation description based on two-body ag- 
gregation eventually breaks down. Our approach is also 
inapplicable to fractal-shape islands, a situation which 
can arise at low temperature JlSj] . Finally, the rate equa- 
tion description should fail below an upper critical dimen- 
sion d c which is given by d c = 2/(1 — /j — A), with 
A the homogeneity degree of the reaction kernel. For our 
point-island model with kernel given by Eq. (Q), A = — // 
and hence d c = 2. Thus the relevant two-dimensional 
case corresponds to the critical dimension, and logarith- 
mic corrections to mean-field predictions, in addition to 
the previously-discussed logarithmic corrections to the 
reaction rate, can be anticipated. 



III. GROWTH WITH FLUX 



Consider a point-island system with flux F of 
monomers onto the substrate. The rate equations for 
the concentrations c& of islands of mass k in the presence 
of this steady monomer flux are 

dc 1 °° 

-j£ = y K H C i C J ~ C kYl K kjCj + F S kl - (2) 



These rate equations represent a mean-field approxima- 
tion in which spatial fluctuations are neglected, and also 
a low-coverage approximation, since only binary interac- 
tions are treated. 
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A. Steady State Regime 

Let us first consider the behavior in the steady state 
regime. To analyze this system of rate equations when 
the time derivative term is set to zero, it proves con- 
venient to use the generating function technique. We 
introduce two generating functions 



DC 



oo 



fe=i fe=i 



(3) 



Multiplying Eq. (El) by z , and summing over all k, gives 



C^(z)C(z) - C^(z)N - C{z)N^ + Fz = 0. 



(4) 



Here N = C(z = 1) = c fc is the total island density 
andA^=C M (z = l) = £fc-"c fc . 

We now assume a power law asymptotic behavior for 
the steady state concentration, 



C 



(5) 



as k — > oo. For this power law to hold for all k, we re- 
quire t > 1, so that k~ T converges; this leads to the 
condition /i < 1 for the mobility exponent as shown be- 
low. From basic Tauberian theorems , the asymptotic 
form for in Eq. (|5|) induces the following power-law 
singularities in the generating functions as z — » 1 

C(z) = N + CT(1 - r)(l - zy- 1 + ..., 
C„(z) = N„ + CT(1 -T- M )(l - z) T+ »- 1 + ..., (6) 

where T is the gamma function. The leading constant 
factor in each line is finite and coincides with the defini- 
tion given in Eq. (^|) if the exponent of the second term is 
positive. Otherwise, the constant factor vanishes and the 
generating function has a power-law divergence as z — > 1. 
Substituting these expansions into Eq. (Q) and matching 
the leading behavior in (1 — z) leads to the decay expo- 
nent t = (3 — A 1 )/2. The condition for a steady state 
to occur, t > 1, thus imposes an upper bound on the 
mobility exponent, /i < 1. From matching the leading 
behavior in (1 — z), the constant C may also be deter- 
mined, from which the island mass distribution in the 
steady-state regime < fi < 1 is 



/ F{1 - n 2 ) cos(tt^/2) , -(3-^/2 
V 4tt 



(7) 



In the special case of constant reaction rate, /i = 0, 
one can find a complete time-dependent solution. In this 
case, the rate equations are 



~dt 



<kCj - 2c k N + F5 kl 



(8) 



i-\-j—k 



By summing these equations over k, one obtains the as- 
sociated equations for the density, 



dN 
~dt 



-AT + F, 



and for the generating function, 
dC(z,t) 



dt 



= C{z, t) - 2C(z, t)N(t) + Fz. 



(9) 



(10) 



Solving these equations for initially clean surface, Cfc(0) 
0, we obtain the time-dependent island density 



N(t) = VF tanh (t^F) , (H) 
and the generating function 
C{z, t) = N{t) - v/F(l - z) tanh [t ^F{1 - zfj . (12) 

Expanding the generating function C(z,t) in a series 
in z, gives the concentrations Cfc(i). Relatively simple 
results are obtained in the long time limit 



l~F T(k 
Ck{t) V 4tt T(k + 1) 



I) 



(13) 



This exact solution of Eq. (|13j) agrees with the asymptotic 
solution of Eq. (]?]) . Note that islands of sizes comparable 
with (1 — z)~ l or smaller, give dominant contribution to 
the generating function, Cq(z, t) — J2k>i zk °k{t), while 
bigger islands provide an asymptotically negligible cor- 
rection. Eq. ( |l2| ) also shows that the crossover from a 
time-dependent to saturated behavior takes place when 
ty/1 — z ~ 1. Together with k ~ (1 — z)^ 1 this im- 
plies that the steady state is established for small islands 
whose size is in the range, k <C t 2 . For general mobility 
exponent, one physically expects a similar crossover to 
steady behavior when k ~ t£, with mass cutoff exponent 
C dependent on ji. To determine £, we use the physical 
condition that the total mass on the substrate kc k (t) 
equal Ft together with the steady state asymptotics of 
Eq. (@). This gives 



(14) 



fe=i 



fc=i 



that is, C = 2/(/x+l). 



B. Continuous Island Evolution 

For sufficiently large mobility exponent, /i > 1, contin- 
uous evolution can be anticipated. Indeed, in the previ- 
ous subsection we showed that assuming a steady state 
leads to fj, < 1. Additionally, continuous evolution is 
known to occur when ji = oo, that is, in the extreme 
case of mobile monomers and immobile islands. This 
case has been treated analytically [||-f§ and power-law 
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growth in the number of islands, N(t) ~ (SF 2 ^) 1 / 3 , was 
found. Moreover, in the scaling region 



oo, t - 



00, 



with x 



2k 



<1, (15) 



(Ft 2 ) 1/3 

the island mass distribution approaches the scaling form 

1/3 



(1-x) 



-1/2 



when x < 1, (16) 



while for x > 1 the island concentrations are negligible. 

We therefore expect that for all 1 < [i < oo, the island 
mass distribution does not reach a steady state. The 
marginal case of p — 1 is interesting since it corresponds 
to experimentally relevant case of terrace diffusion [ jl3| . 
As discussed in the previous section, such values of the 
mobility exponent p naturally appear for different micro- 
scopic mass-transport mechanisms. 

Our primary results are that when p is strictly greater 
than unity but still finite, 



N(t) 



sin(7r//x) 



In (tyfF 



M/2 



(17) 



while the concentration of islands of mass k decays in 
time as 



c k (t)~VF(k\) 



In ( tVF 



-ju(2fc-l)/2 



(18) 



for sufficiently small islands, k <C ln(tvF). It is remark- 
able that these logarithmic dependences, a feature which 
generally signals marginal behavior, occurs in the entire 
regime 1 < p < oo. 

Our approach that leads to Eqs. ( Jl7| ) and ([ll]) is based 
on the physical picture that the system is slowly evolving 
because the growth of islands by aggregation is substan- 
tially counterbalanced by the input of monomers. These 
competing effects lead to nearly time-independent island 
concentrations over an "inner" size range, while more 
strongly-time-dependent behavior occurs in an "outer" 
range which extends to the largest islands. In the inner 
region, the picture of near balance between aggregation 
and input motivates the use of the quasi-static approxi- 
mation, where the time derivative in Eq. (|^) is initially 
neglected, from which the island concentrations as a func- 
tion of the total concentration of islands can be obtained. 
The dependence of the island concentrations on time is 
then determined by the condition that the total mass in 
the system is proportional to t. The validity of this ap- 
proach may be verified a posteriori, where the logarith- 
mic dependences in Eqs. (0) and imply that the 
temporal derivatives in the Smoluchowski rate equations 
are asymptotically negligible. 

Within this quasi-static framework, Eqs. (^) become 

= 1 - ci (N + Ay , 

1 



By summing Eqs. ( |l9| ) over all k, the total island density 
in the quasi-static limit obeys 



= 1- NN U 



(20) 



Eqs. (H) and © have been non-dimensionalized by the 
scale transformation c k — > y/F c% . We also scale the time 
variable by t — > tj \[F, so that the mass density obeys 



6(t)=J2kc k (t) 



(21) 



fc=i 



Eq. ( |20| ) immediately gives = N 1 , and then from 
the first of Eqs. @, c\ ~ 1/N. The remainder of 
Eqs. (|l^) may then be solved recursively. By writing 
the first few of these equations, it is evident that the 
dominant contribution to c k is the term in the quadratic 
product which is proportional to c\c k -\. Keeping only 
this contribution, the resulting recursion may be solved 
straightforwardly to yield 



Ck ^N 



n i + jM/AT2 

vll" IK 



ff(i+r") 

3=1 



N 



(22) 



3=2 3=1 



Since the factors B 3 < 1 for < N 2 , while Bj — > 1 
for ^> A^ 2 , this implies that c k is a rapidly decreasing 
function of k for k <C N 2 ^ = k and then converges to a 
finite value p for k > n. 

To compute c k , first note that for /i > 1 the product 
Ylj bj converges, so that it may treated as constant. We 
then write the second product as the exponential of a 
sum and take the continuum limit. This leads to 



Cfe 



1 




N 


exp 


1 




N 


exp 



\l+j»/N 2 



-N 2 ^ 1 / ln(l + w->*) dw 
Jo 



(23) 



where w — j/N 2 ^ and x = k/N 2 ^. This form has 
two slightly different asymptotic behaviors depending on 
whether /i is strictly greater than or equal to 1. For 
fi > 1, the monotonically increasing integral in Eq. ( [23] ) 
converges as x — * oo. Thus c k decreases as a function 
of k until a threshold value k ~ N 2 ^, beyond which c k 
remains constant with a value p determined by taking 
the upper limit of the integral as infinite. Hence 



('k 



A7 exp 



(24) 







2 ? 

i-\-j—k 



+ Ci Cj - c fc (fc-^A^ + N^) . (19) with An = / °° ln(l + w^) dw = tt/ sin(7r//x). Within the 

quasi-static approximation, this constancy in c& should 
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persist over the range, n < k <C K. Here K defines an 
upper limit for the "inner" regime where the quasi-static 
approximation remains valid. Physically, this limit corre- 
sponds to the size range where islands are only beginning 
to form. 

Because of the importance of the temporal decay of 
the densities in the outer size range, k 3> k, the quasi- 
static approximation is inadequate and an alternative is 
needed. We therefore separately account for the distribu- 
tion of these "raw" (evolving) islands and then perform 
a matched asymptotic expansion p2| , to join the inner 
(k <C K) quasi-static solution of "ripe" islands to the 
outer (k 3> k) solution of raw islands in the overlap re- 
gion k <C k <C K. In the outer region, the asymptotically 
dominant terms in the rate equations are P3[ 



dc k . . 
-— = ci (c fc _i - CfcJ . 
at 



(25) 



Since raw islands are large, we employ the continuum 
limit of Eq. pi), 



d did 2 

ar + 9fc~2ap ) Cfe(T) = ' 



where T = J c \dt' c x {t'). If wc neglect the diffusive term, 
the solution of the resulting wave equation is Cfc(T) = 
f(T — fc), with / an arbitrary function. Matching with 
the inner solution determines / and gives 



c k {T)=p{T-k), 



(27) 



which thus provides the time dependence of the is- 
land concentrations once T has been determined. Also, 
Eq. ( p7j ) provides the estimate K ss T for the cutoff size 
beyond which essentially no islands exist. 

To close the solution, we now determine the time de- 
pendence of the island density N(t). To accomplish this, 
we use the sum rules for the mass and island densities, 
^] kck = M = t and ^2 Cfc = N. In the continuum limit 
these become 



t 



T 

du (T — u) p(u) , N = I dup(u) 
Jo 



(28) 



The dominant contribution to these integrals arises from 
large (raw) islands in the outer region. The second rela- 
tion also gives ^ = p(T). Combining this with Eq. ( |2^ ) 
yields 



P(T) 



(InT) 



ft/2-1 



T 



(29) 



The raw island mass distribution is then given by 
Eqs. @ and <M) as, 



c fe (T) 



[ln(T-fc)] 



H/2-1 



T-k 



(30) 



out by the diffusive term in Eq. (|2^) . Thus instead of the 
singularity at k = K , the density of raw islands reaches 
a peak value of the order of t~ x l 2 and then rapidly de- 
creases for larger k. 

Using Eq. (|29j), we can now estimate the left integral in 
Eq. ( p8| ) and find that the ratio of the first to the second 
term scales as InT. Keeping only the dominant contribu- 
tion and using the right integral of ( |2§| ) gives T ~ t/N. 
Thus the cutoff size is given by 



t t 
{t >- N~ (lnt)^/ 2 ' 



(31) 



Parenthetically, the basic relation between the real and 
modified times, T = J* dt'ci(t'), together with c\ = 1/N 
also leads to T ~ t/N and thus to Eq. (31). 
In real time the island density becomes 



N(t) =s A- 1 In 



lnt 



l-3p/2' 



H/2 



(32) 



(26) where = it/ s'm(ir/p). Upon neglecting the logarith- 
mic temporal factor inside the logarithm, our basic result 
quoted in Eq. (frh is recovered. 



C. The Case fi = 1 

In the borderline case p, = 1, subtler nested logarithmic 
behavior arises, as reflected by a singularity in Eq. jl^ ) 
that appears upon formal continuation to p — > 1. The 
analysis closely parallels the case p > 1 and we just out- 
line the main results. 

In the inner region, the general approach of the pre- 
vious subsection applies up to Eq. (|22j). However, when 

" l + .r 1 ) inEq. M) 



/' 



1, the product flt=i 



rfc-l, 



,-, , ii: ri • ./ 

now equals k, i.e., it diverges. Second, the term C2Ck-2 in 
addition to c\Ck-\ contributes to the asymptotic behav- 
ior. Due to this latter attribute, the recursion relation 
for Cfc becomes 



Cfc 1 
k 



k/N 2 



k/N 2 



k-l 



Cfc-2 



1 

2 N 2 



TrI- (33) 



We seek a solution for Ck in the form of Eq. (|22j) . Thus 
we write 



a,~C k - 



i=2 



j/N 2 



■j/N 2 



(34) 



where the factor Cb accounts for the additional term in 



Eq. (33). Substituting into Eq. (33) gives the recursion 
formula for this correction factor Cfc 



Notice that the singularity when k — K = T '. However in 
the mass range T — k ~ \T, this singularity is smoothed 



Cfc — Cfc_i + Cfc_2 



1 



k - 1 



1 

N 2 



(35) 
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These coefficients are slowly varying in k when k ^> 1 and 
we may treat k as continuous in this asymptotic regime. 
Eq. (J3^) then becomes a differential equation whose so- 
lution is Cfe ~ ke x (with x = k/N 2 ). Consequently, 



N 



exp 



x-N 1 



lnfl 



L ) dw 



(36) 



Thus for /i = 1, the concentration Cfe decreases rapidly in 
k for k <C ./V 2 and reaches a minimum at fc = iV 4 whose 
value is 



exp [-N 2 In TV 2 ] 



(37) 



In the outer region, the general results of the previous 
subsection, Ens. ( p5| ) through (|28|), are still valid. Thus 



combining (|37|) with jjrf = p(T) we solve for the total 
island concentration to obtain 

N(t) 



hit 



rn(lni) 



Eq. ([34|) together with Ck ~ fc implies that the concen- 
tration of islands for mass k <C A^ 2 is 



Ck(t) 



(fc + 1)! 



(fc + 1)! 



ln(lnt) 



-i fc-l/2 



lni 



(39) 



Then the island mass density is approximately constant, 
Cfc(T) ~ p(T), where 

p(T) - T-^ln^-^lnQnT)]- 1 . (40) 

This holds over the range k < < K with K(t) ~ 
t ^/ln(lnf)/lnt. Finally, the raw island density given by 
Eqs. ( p7| ) and ( |40| ) holds up to a mass cutoff 

D. Maximum Island Density 

To apply our results to real submonolayer epitaxial 
systems, note that in the submonolayer regime, the cov- 
erage must be small, that is, M = Ft -C 1. On the 
other hand, the asymptotic predictions of our theory ap- 
ply when ty/F 3> 1. Consequently, our results should be 
valid in the time range F~ x l 2 <C t <C F^ 1 . Since the 
(dimcnsionless) flux F is typically small in epitaxy ex- 
periments, the time range over which our theory should 
apply is correspondingly large. A commonly employed 
connection between theory and experimental results is to 
determine the maximum island density at the end of the 
submonolayer regime, £ m ax ~ F . The conventionally- 
quoted result is that this maximum density scales as a 
power of the flux p|,p4|] 



F x 



(41) 



with the flux exponent typically in the range ^ < X < 
~ ^HlUm. While most theoretical studies predict 



X < 1/2 p| Ji^ , p7| , p8| , (e.g., for point immobile islands 
X = 1/3), the value \ = 1/2 has been observed experi- 
mentally, see Ernst et al. p6[ . Our analysis predicts that 
the maximum island density attains the value 



iVn 



F 1 / 2 [\n{l/F)] 



M/2 



(42) 



and that \ = 1/2 is generic and applies to any model 
where island diffusion leads to continuous evolution. Pre- 
vious theoretical work has focused on somewhat patho- 
logical models where only a few (the smallest) island 
species could diffuse. In such models, the exponent x 
is sensitive to the cutoff size between mobile and immo- 
bile islands. Fortunately, the generic situation is simpler 
and the asymptotic flux dependence of N maDt is univer- 
sal. However, there do exist logarithmic factors which 
are non-universal, as they depend on the mobility expo- 
nent. For example, the logarithmic factor in Eq. ( |42| ) 
gives an effective exponent, when plotting N max versus 



(38) F, XefiF 



1 



ln(ln(l/F)) 
I 1 ln(l/F) 



which is smaller than 1/2. 



Indeed, for a fixed flux, the effective flux exponent is a 
decreasing function of fi. Such a feature should be taken 
into account in the interpretation of experimental and 
numerical data. 



IV. GROWTH WITHOUT FLUX 

We now consider epitaxial growth for point islands 
with a non-zero initial monomer density, reaction rate 
Kij = i~^ L + and no subsequent monomer flux. Such 
a system has been extensively investigated within the 
framework of irreversible aggregation |29| -^3|], as well as 
in theoretical studies of point-island aggregation with im- 
mobile islands ||,[34j . More realistic examples of the lat- 
ter system have also been investigated numerically and 
experimentally in recent studies (see p^j35[| and refer- 
ences therein). 

We give here a simple argument for the asymptotic 
form of the island size distribution. This argument is 
based on first solving for the island and monomer densi- 
ties, N(t) and Ci(i), respectively, and then using scaling 
]l7| to infer the asymptotics of the distribution. Accord- 
ing to the scaling ansatz, the asymptotic island-mass dis- 
tribution should have the following form 



c k (t) ~ N 2 G(x), x = kN, 



(43) 



for finite x. The constraints JdxG(x) — 1 and 
J dx x G(x) = 9 automatically enforce the conditions 
J2 c k(t) = N(t) and mass conservation J2^°k(t) = @- 
Starting with the exact rate equation for N(t), 



dN 
~dt 



1 oo oo 
i=l j=l 



' Ci Cj 



(44) 



we substitute the scaling ansatz Eq. ( p3[ ) to obtain N 
_]\j-2+fi w j lose solution is 
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N(t) ~ f-VCi+M). 



(45) 



Consider now the rate equation for the monomer den- 
sity, ci = -Kci(N + N,,). As shown by Eq. @, we 
may neglect the second term and integrate the resulting 
equation to yield 



ci (t) ~ exp 



_ t n/{l+n) 



(46) 



Together with Eqs. (|43|) and (^5j), (|6| predicts that the 
asymptotic behavior for the scaling function in the small- 
s' limit is exp(— l/x M ) . On the other hand, for the entire 



class of reaction kernels of the form = i 



3~ 



the scaling function G(x) decays exponentially |17J in x 
for large x. Combining these results gives the asymp- 
totic forms of the island size distribution in the absence 
of monomer flux 



G(x) 



-l/x" 



x[ 0, 
x ] oo. 



(47) 



To determine the validity of these mean-field predic- 
tions, let us consider the kinetics of this process in gen- 
eral spatial dimension. To this end, we determine the 
reactivity of an arbitrary cluster by a dimensional argu- 
ment. Consider an island of radius R which diffuses with 
the diffusion coefficient D in d-dimensional space. Dur- 
ing a time interval t this particle traces out the so-called 
Wiener sausage whose volume is |3(J 



v d (t) 



Dt [\n(Dt/R 2 )]' 
DtR d ~ 2 , 



d=2, 
d>2. 



(48) 



For epitaxial growth with no flux, R corresponds to a 
(growing) average island radius, and D to the diffusiv- 
ity of a mass-fc island, D ~ fc~ M ~ R~ dfl . Clearly, all 
monomers initially within a Wiener sausage will aggre- 
gate into a typical island by time t. Consequently, the 
initial coverage of the substrate 9 is given by R d ~ 6Vd(t). 
This, together with Eq. gives R 2 ~ Dt in all dimen- 
sions d>2. Finally, combining R 2 ~ Dt and D ~ R~ dfJ - 
leads to 



(49) 



In two dimensions, the rate equation predictions 
(Eq. J45|)) and the above heuristic argument agree. Thus 
in the absence of a monomer flux, and under the assump- 
tion of compact islands, the two sources of logarithmic 
corrections - a logarithmic dependence of the reaction 
rate on island radius and the fact that the system is 
at the critical dimension - evidently cancel each other. 
Therefore the Smoluchowski rate equations appear to be 
asymptotically correct in two dimensions. 



V. NUMERICAL RESULTS 

We have performed detailed Monte Carlo simulations 
for island growth for two models of epitaxial growth. We 



first consider point islands, in which single-site islands 
hop with equal probability to any lattice site (diffusion 
on a complete graph) and aggregate whenever two islands 
occupy the same site. This corresponds to the mean- 
field limit and thus provides a direct test of some of the 
delicate approximations made in our theoretical analy- 
sis. The second model is a more faithful description in 
which circular island "droplets" with radius proportional 
to the square-root of the island mass diffuse to nearest- 
neighbor sites on a two-dimensional lattice. Whenever 
there is overlap of two islands, they immediately coalesce 
to a single island which is centered at the initial posi- 
tion of the larger island. After each coalescence, a test is 
made to determine if additional overlaps have been cre- 
ated. All such higher-order coalescences are performed 
until all overlaps are resolved. We treat here only the 
case with external flux, since growth without flux has 
already been investigated numerically (see e.g., p3[ and 
references therein). 

The point-island model is relatively easy to imple- 
ment. Additionally, this model has a technical advan- 
tage over growing droplets in that the submonolayer de- 
scription applies over a longer time range. For point 
islands, the submonolayer regime is defined by the cri- 
terion N <C 1, which is considerably less stringent than 
the criterion Ft <C 1 appropriate for the droplet-island 
model. Therefore our theoretical results may be com- 
pared with the point-island simulations in the time range 
F~ x l 2 <f< exp (—i* 1-1 / M ) and with droplet-island sim- 
ulations in the time range F -1 / 2 <C i <C F^ 1 . 



A. Point Islands 

Due to the point-like nature of islands and the 
equiprobable hopping to any site of the system, the sim- 
ulations correspond directly to the Smoluchowski rate 
equations. Simulations of point islands were performed 
on a graph of L sites with the following time evolution. 
At any stage, the number of deposition events per unit 
time is FL 2 , while the number of aggregation events per 
unit time is L 2 CC . Here L occ is the number of occupied 
sites in the system. In a microscopic event, deposition 
is chosen with probability rd ep = FL 2 /(FL 2 + L 2 CC ) and 
aggregation is chosen with probability r agg = 1 — rd cp - 
In a deposition event, a monomer is deposited onto a 
randomly chosen vacant site. We checked that restrict- 
ing the deposition onto unoccupied sites only does not 
alter our results. In an aggregation event, two sites were 
chosen randomly from the list of occupied sites. If these 
sites contain islands of mass i and j, aggregation occurs 
with probability (i~ M + j _A1 ) /2. Time is then increased 
by At = 1/(FL 2 + L 2 CC ), and the procedure is repeated. 

To test our algorithm, we considered first the extreme 
case of immobile islands, \x — 0, where a complete time- 
dependent analytical solution is available (|l^). Simula- 
tions were performed for an initially empty system with L 
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between 10 3 and 10 4 . We found excellent agreement be- 
tween numerical and theoretical predictions, both for the 
steady state and transient characteristics. More generally 
for \i < 1, our simulations showed that the system reaches 
a steady state with the observed steady state character- 
istics in good agreement with the theoretical prediction 
ofEq. flf. 



island density (Fig. 2) are consistent with N(t) growing 
as a power of Int, as predicted by Eq. (0), but with 
somewhat smaller exponent than /i/2. We consider our 
data sufficient to exclude a power-law time dependence 
of the island density. However, impractically long simu- 
lation would be needed to determine the exponent of the 
logarithmic factor in Eq. (|l7|). 



o 

-2 




FIG. 1. Semilogarithmic plot of lnc k (t) vs k at t « 22000 
for n = 1.5 based on simulations of point islands. The data 
is based on 5000 realizations of an initially empty system of 
2000 sites with F = 0.05. Notice the existence of an inner 
scale k < k w 50 (inset) and an outer scale for k « K m 4000. 



4 




log(tF ) 

FIG. 2. Plot of N(t)/VF vs. log (frv/F) for fi = 1.2 (O), 
fi = 1.4 (★), and fi = 2.0 (□). The simulation parameters are 
the same as in Fig. 1. 

For /j, > 1, the total number of islands is found to in- 
crease indefinitely and the island mass distribution sub- 
stantively agrees with our theoretical predictions (see 
Fig. 1). Indeed, the island mass distribution Ck(t) sharply 
decreases for small mass (ripe islands), and then increases 
in a mass range which corresponds to raw islands. In the 
proximity of k rs K, there is a peak in cj, as predicted 
by our description based on matched asymptotic expan- 
sions. Our data for the time dependence of the total 



B. Circular Island Droplets 

We now consider simulations of compact growing cir- 
cular islands which we term as droplets. We additionally 
assume that islands are always centered on sites of the 
square lattice. In the simulations, we consider a system 
of size L 2 to which we add monomer droplets of radius 
ro = 0.495 to guarantee that adjacent monomers do not 
overlap. Monomers, however, can overlap with bigger 
droplets. When two droplets of radii r\ and T2 overlap 
they coalesce to form a droplet of radius \Jr\ + r\ which 
is located at the center of the larger of the two initial 
coalescing droplets. As the coverage increases, multiple 
coalescences become increasingly probable and we there- 
fore check for new island overlaps after each coalescence 
event and perform additional coalescences and continue 
to resolve all additional overlaps, if needed. 

In the time evolution, a microscopic process, either 
deposition onto any system site or diffusion to a nearest- 
neighbor site is chosen with respective probabilities 
Pdop = F/(F + N) and p dif = 1 - p dcp . We then test 
for and perform all possible subsequent coalescences af- 
ter each event. After completion of this microscopic event 
the time is incremented by At = (FL 2 + NL 2 ) -1 and the 
process is repeated. 

2.0 1 




FIG. 3. Plot of N{t)/V~F versus tVF for fx = 1.2 and 
F — 0.001 for point (□) and growing circular islands (*). 

Fig. 3 compares the time evolution of point and grow- 
ing circular islands. For t < 1/yF, the plots of N(t) 
versus t for the two processes coincide, thus indicating 
that the point-island model provides an excellent early- 
time description for the more realistic model of growing 
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droplets. However, when t > 1/yF, multi-body aggre- 
gation starts to become important and the density of 
droplets decreases with time while the density of point 
islands continues to grow. Thus the point-island model 
is a suitable starting point to interpret simulational and 
experimental data of epitaxial growth |38[ . 

A visualization of the aggregation of circular droplets 
at a relatively late stage is shown in Fig. 4. Here, many- 
body effects become significant and a coalescence of two 
large droplets can lead to a large avalanche of coalescence 
events. 



c 



Q, n . 



OPS 
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FIG. 4. Central portion of a typical droplet configuration 
for fi = 1.2 and F = 0.001 for a system of size 500 x 500. 



We have also performed simulations of two- 
dimensional fractal islands. Specifically, when two 
monomers occupy neighboring sites they stick irre- 
versibly to form a dimer. This process continues and 
leads to the formation of fractal islands which hop as a 
rigid unit. Simulation of this process is simpler than in 
the growing droplet model because coalescences subse- 
quent to the primary event cannot occur. This model has 
been studied earlier [|l4| and our respective results agree. 
For example, for the case [i = 1.2 we obtain x ~ 0.43 if 
we fit the data for -/V max versus F to power-law behavior. 
Our interpretation for this exponent value differs, how- 
ever, from that of Ref. |L4|]. Indeed, they suggest that 
there might be a deep connection to a point-island model 
where both monomers and dimers diffuse and larger is- 
lands are immobile, since the exponent in that model, 
X = 2/5 |27l], is close to their simulation results. We 
believe that there is no connection between these two 
models and the exponent observed by simulation is just 
an effective value whose asymptotic valu e is governed 
by the logarithmic correction in Eq. (B3). Our theory 
predicts universal flux exponent x — 1/2) and the dif- 
ference with observed effective flux exponent is due to 
the logarithmic correction. When this feature is taken 
into account, a value for x very close to our prediction is 
found. 



VI. SUMMARY AND DISCUSSION 



In Fig. 5 we plot the dependence of the maximum is- 
land density as a function of the flux. Using the power- 
law form of Eq. ( f4l| ) leads to a best-fit value x ~ 0.43, 
while taking into account the logarithmic correction given 
by Eq. ([|2|) gives x ~ 0.50, in excellent agreement with 
our theory. However, using the additional flux renormal- 
ization from our more accurate treatment of the reaction 
rate (Appendix B) to fit the data leads to the somewhat 
different exponent estimate of x ~ 0.53. 
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FIG. 5. Plot of log(Ar max /[log(l/f)] (1+fl)/2 ) versus log(F) 
for n — 1.2. The slope of the line is 0.53. 



We have investigated the kinetics of submonolayer epi- 
taxial growth within a simple model which incorporates 
basic physical features of epitaxial growth - deposition, 
island diffusion, and aggregation. We have shown that 
our model displays universal kinetics - up to logarithms 
- as the governing exponents associated with the time 
and flux dependence of observables are detail indepen- 
dent. This is in contrast to the behavior exhibited by 
models where islands below a cutoff size are mobile while 
larger islands are immobile. In this case, characteristic 
exponents depend on this mass threshold. 

We analyzed in detail the situation where the effective 
diffusivity Dk for islands of mass k is cx fc _A \ Such 
a diffusivity arises, e.g., in periphery diffusion, where an 
adatom on the edge of an island detaches, hops to a neigh- 
boring site on the edge, and then re-attaches to the is- 
land. This mechanism also causes islands to be compact. 
A Smoluchowski approach shows that the reaction rate 
between two islands of mass i and j is Kij cx + 
multiplied by a factor which depends logarithmically on 
their radii in two dimensions. This weak dependence im- 
plies that a model with point-like diffusing islands should 
be quantitatively accurate when applied to the submono- 
layer regime. The net effect of the logarithmic factor 
is the flux renormalization, F — > Fhi(l/F), as demon- 
strated in Appendix B. 

For mobility exponent < fi < 1, there is a steady 
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state, with the concentration of islands of mass k, Ck oc 
fe~( 3 ~w/ 2 . For (i > 1, logarithmic time dependence arises 
in which the total island density N(t) oc (ln£) M / 2 . In this 
regime, the island distribution exhibits a rich mass de- 
pendence in which there is: (i) a precipitous decay in a 
"boundary layer" k <C k (with k ~ lnt), (ii) a grad- 
ual growth in the main part of the mass distribution 
k < k < K (with K ~ t(lnfT^ 2 ), and (hi) an in- 
ternal layer \k — K\ ~ \ft where the density of islands 
reaches a peak and then sharply vanishes. The entire 
regime 1 < /i < oo exhibits this same behavior up to 
logarithmic corrections, while the transition between the 
steady and evolving regimes at /i — 1 is characterized by 
nested logarithmic behavior. 

It is noteworthy that our theoretical approach can 
be applied to epitaxial systems with arbitrary mass- 
dependent diffusivity D k which decays faster than inverse 
mass. For this general situation, the analog of Eq. (|2S" 
is 



k 

"•Jin 



_ + N 2 Di ' 

2 J 



(50) 



For example, for a diffusivity which decays exponentially 
in island mass, D k = e _a ( fc_1 ), the case investigated nu- 
merically in Ref. |39], our theory predicts 



N(t) 



Fexp 



| In (tVF 



(51) 



and K ~ t/N, leading again to the same universal value 
of the mass cutoff exponent £(a) = 1. Eq. (51) exhibits 
unusual time dependence - faster than any power of log- 
arithm and slower than any power law - and may be 
difficult to observe numerically. The maximum island 
density is 



Fexp 



ln(l/F) 



(52) 



so again \ = 1/2- Numerically, the exponent x( a ) a P~ 
pears to decrease as a increases [§9|. Our analysis sug- 
gests that in the asymptotic regime x( a ) = 1/2 for 
all < a < oo. However, fitting the functional form 
in Eq. (p^) to a single power-law in F gives Xcfi — 
\ - y/a/[4 ln(l/F)]. Therefore even for small flux the 
effective exponent may be considerably smaller than 1/2. 
Also, Xeff(o) is a decreasing function of a, in agreement 
with the observations from the simulation | p9[ . 
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APPENDIX A: EPITAXIAL GROWTH ON ID 
SUBSTRATES 

We now extend our results to the case of a one- 
dimensional substrate. Since the upper critical dimen- 
sion d c = 2 po|| , the mean-field approximation does not 
apply in one dimension. In the absence of a theoretical 
framework to systematically treat the case d < d Cl we 
give a heuristic treatment. We will derive results for the 
specific cases of d = 1 and d = 2; comparison between 
the latter and the rate equation results provide a check 
of our approach. 

Consider first the simpler case of systems that ap- 
proach a steady state. We present an argument based 
on the volume swept out by a Wiener sausage, as in the 
case of systems without flux (Sec. IV). To mimic the ef- 
fect of the flux, we suppose that there is no flux but that 
all islands have initial mass which equals t. At time t, all 
islands within a reaction volume (Dt) d / 2 have coalesced 
into a single island. Ignoring logarithmic correction in 
two dimensions, this gives the following estimate for the 
average island mass M, 



M(t) - t x 



Dt, 
VDi, 



(Al) 



Combining (Al) with D ~ Af _M and using the fact that 
the average island mass M scales as the mass cutoff we 
find 



K{t) 



t 2/(l+M) > d = 2, 
<3 /(2+ M ) d = h 



(A2) 



In the steady state regime Ck ~ k T . Consequently the 
sum rule t = kck ~ K 2 ~ T ~ t^ 2 ~ r \ im plies the rela- 
tion t = 2 — This, together with (A2), gives 



(3-m)/2, d = 2, 
(4-m)/3, d = l. 



(A3) 



This argument reproduces the correct values of the cut- 
off and the decay exponents when d = 2, and we an- 
ticipate that the one-dimensional results are also exact. 
Indeed, an exact solution Q of one-dimensional aggre- 
gation with monomer input and mass-independent island 
diffusivities gives r = 4/3. Finally, note that the expo- 
nents C and r attain the critical value £ = r = 1 when the 
mobility exponent /i = 1 in both one and two dimensions. 

When fi > 1, continuous evolution occurs. We now 
argue that the densities asymptotically evolve according 
to generalized rate equations. For monomers, we write 



dci 
~~dt 



F 



At' 



(A4) 



Here At is the collision time for a monomer to encounter 
an island. During this time interval, a monomer visits 
yAt different sites in one dimension, so that the colli- 
sion time is determined by NyAt 1. Consequently, 
the rate equation for the monomer density becomes 
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at 



(A5) 



Continuing this line of reasoning we obtain the rate equa- 
tions for the densities Ck(t) which differ from Eqs. (||) by 
a factor TV in each reaction term. 

We now analyze these equations by the same quasi- 
static framework as in Sec. II. Thus we need to solve 

= F-iV 2 ci, 

° = 5 51 (i-^+j-^CiCj-Ckik-^N + N^). (A6) 

i+j=k 

Repeating the steps of our previous derivations we ob- 
tain, e.g., for the density of islands, 



N(t) ~ F 1 / 3 



sin(7r//u) 



In 



fi/2 



(A7) 



and for the behavior of the density of relatively small 
islands 



c fc (t) ~ F 



1/3 



(A8) 



This in the continuously evolving regime, the time depen- 
dence remains primarily unaffected by the dimensional- 
ity of the substrate. However, the flux dependence does 
change with d, and we find N max ~ F 1 / 3 Jhi(l/F)]^/ 2 
(compared to the F 1 / 2 dependence in Eq. (|42|)). Overall, 
universal behavior arises in continuous evolution which is 
only slightly affected by model details, substrate dimen- 
sionality, etc.. 



APPENDIX B: REACTION RATE IN TWO 
DIMENSIONS 

We now account for the logarithmic corrections to the 
reaction rate (Eq. (jlj)) that appear in two dimensions. 
We first demonstrate that these logarithmic corrections 
can be accounted for within modified rate equations. Let 
us first consider a simpler model where point islands dif- 
fuse at the same mass-independent rate. Then the total 
island density c(t) obeys c = —c/At, where At is the 
time between successive collisions. A collision is expected 
when the island visits 1/c distinct sites. Since the num- 
ber of distinct sites visited by a random walk in time t 
grows as Dt/\w(Dt) in two dimensions |3(|, the collision 
time follows from the condition DAt/\n(DAt) ~ 1/c. 
The resulting expression for At leads to the rate equa- 
tion c = ~Dc 2 j ln(l/c). Similarly, for an island of radius 
R we obtain c = -£>c 2 /ln(l/ci? 2 ). 

For growing droplets with mass dependent diffusivity, 
these logarithmic factors imply that the reaction rate 
Kij = Di + Dj should be replaced by 



K, 



A + D 4 



]n[N-i(Ri + Rj)- 2 } 



(Bl) 



In the low coverage limit, the average separation between 
neighboring islands A -1 / 2 is much larger than the aver- 
age island size. Keeping only this dominant factor in- 
side the logarithm gives the asymptotic form of the reac- 



tion rate ~ (Di + Dj)/ln(l/N). Moreover, we can 
replace the total island density by y/F inside the log- 
arithm. This is obvious when ji < 1, since in this case 
the island density indeed approaches a steady state value 
Noo ~ v~F. For (j, > 1, the island density grows according 
to N ~ \fF(\nt)^/ 2 . However, the time-dependent fac- 
tor is clearly sub-dominant as it is at most logarithmic in 
the flux, (lnimax)^ 2 = (ln(l/F))^/ 2 . Hence for arbitrary 
mobility exponent /x, the form ~ (Di + Dj)/ ln(l/F) 
provide an asymptotically correct description for the re- 
action rate. 

Therefore in two dimensions we can continue to use the 
mean-field Smoluchowski equations, with the modifica- 
tion of the reactive term by the factor l/ln(l/F). Upon 
rescaling the densities by ct — > Ck\/ F\n(l/ F), and the 
time variable by t -> txjF^ 1 ln(l/F), we formally map 
Smoluchowski equations for epitaxial growth in two di- 
mensions onto the mean- field Eqs. (^) with F = 1 and 
the reaction rates given by . 

We therefore conclude that we can apply the mean- 
field results to two-dimensional substrates upon making 
the flux renormalization 



F -> Fln(l/F). 



(B2) 



This renormalization does not alter the basic predictions 
of the Smoluchowski approach; for example, all expo- 
nents remain the same. However, this renormalization 
does alter some logarithmic factors, e.g., Eq. (42) for the 
maximum island density is replaced by 



iV max ~f 1/2 Kl/F)](" +1 W 2 



(B3) 
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